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The Arithmetic Cosine Transform: Exact and Approximate Algorithms 


R. J. Cintra* V. S. Dimitrov^ 


Abstract 

In this paper, we introduce a new class of transform method — the arithmetic cosine transform (ACT). We provide the central 
mathematical properties of the ACT, necessary in designing efficient and accurate implementations of the new transform method. The 
key mathematical tools used in the paper come from analytic number theory, in particular the properties of the Riemann zeta function. 
Additionally, we demonstrate that an exact signal interpolation is achievable for any block-length. Approximate calculations were 
also considered. The numerical examples provided show the potential of the ACT for various digital signal processing applications. 
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1 Introduction 

The arithmetic Fourier transform (AFT) has been emerged in 1988 as a signal processing tool (TJ. In a broad sense, the AFT is a 
number-theoretic algorithm for spectrum evaluation. One of the main features of the AFT is its virtually multiplication-free nature. 
Although initial versions of the AFT procedure could only compute the real part of DFT IIE], additional improvements due to Reed 
et al. expanded the AFT methodology and the DFT could be fully evaluated eu. 

Arithmetic methods lack a larger degree of adoption mainly because data is required to be sampled nonuniformly. In fact, the 
AFT algorithm imposes a sampling scheme that can be related to Farey fractions 0. This can be a sensitive issue since incoming 
signals are usually uniformly sampled. Therefore, in order to obtain the necessary nonuniform samples, literature essentially 
suggests two methods: (i) signal oversampling or (ii) sample interpolation. On the one hand, oversampling can be a major drawback 
because it requires an above Nyquist rate sampling. Thus, this is often discarded as an option 0. On the other hand, to obtain 
nonuniform samples from uniformly sampled data, interpolation methods have been considered. Archived literature lacks an exact 
procedure to interpolate nonuniform samples obtained from arithmetic transform methods. Therefore, in a simplistic way, zero- 
order approximations and linear interpolation methods have been considered 0. Although not furnishing exact computations, these 
crude interpolation methods could attain acceptable approximations, when large block-lengths were considered ffl. However, for 
small block-lengths, the implied interpolation errors could be large enough to totally preclude a meaningful computation. 

Nevertheless, recent existing works have applied the AFT algorithm in a number of ways. In Q, the AFT is considered as an 
alternative to the Goertzel algorithm for the single spectral component evaluation. The AFT could also provide frequency domain 
testing units for built-in self-test routines with reduced hardware overhead 0. Additionally, hardware considerations have been 
directed to enhance the associated nonuniform sampling of the AFT 0- Finally, the AFT has been considered as a tool for DCT 
evaluation ma¬ 
in all above mentioned applications, the AFT procedure as described in ||T][3) was considered. In particular, even when the DCT 
was required as in flOl . it was obtained by means of the AFT |j2). Indeed, the DFT spectrum can be mapped into the DCT spectrum 
at the cost of extra computations. 
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^V. S. Dimitrov is with the Advanced Technology Information Processing Systems (ATIPS) Laboratory, University of Calgary, Calgary, AB, Canada. This work 
was done during second author’s sabbatical leave at Nanyang Technological University, Singapore. Email: vdvsdl03@gmail.com 
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The goal of this paper is twofold. First, we aim to introduce a class of purely arithmetic algorithms for the DCT, termed 
arithmetic cosine transforms (ACT). Such new methods are the main contribution of this work. To accomplish this objective, we 
examine arithmetic averages based on the existing arithmetic Fourier transform. Then, we advance new arithmetic averages specially 
tailored for the DCT computation. 

Second, an exact interpolation method for the ACT is sought. Such exact interpolation could provide solid mathematical grounds 
for arithmetic transform methods, constituting a theoretical advance in the area. In fact, any arithmetic transform method could only 
furnish an exact computation if a precise interpolation method could be devised. In particular, we discuss an interpolation scheme 
that could provide an exact spectrum evaluation, even when short block-lengths are considered. Moreover, an asymptotic analysis 
for the proposed interpolation approach is elaborated and a heuristic approximation algorithm is devised. 

The article unfolds as follows. In Section [5] we examine some of the existing tools employed in the AFT theory and propose 
a new algorithm for the DCT. Section [3 is devoted to the interpolation issues that arise from the introduced ACT. In Section [4] the 
discussed arithmetic transform is formatted in terms of matrix representation. Conclusions and final remarks are given in Section^ 


2 Arithmetic cosine transform 


2.1 Preliminary concepts 

Among the several types of DCTs, we separated the DCT-II, which can be regarded as the most employed form HD- This transfor¬ 
mation relates two /V-dimensional vectors, v and V, according to the following pair of relations 03 


V k 



nk(n + \/2)\ 
N )' 


k = 0,l,...,N- 1, 


( 1 ) 


N -1 


= \/ — Yj a k V kCOS 
^ k=0 


nk(n +1 /2) 
N 


n = 0,1,... ,1V— 1, 


where the coefficients a k , k = 0,1 ,...,N — 1 are given by 


1 /s/2, iffc = 0, 

<4 = < 

I 1, otherwise. 

Hereafter, we refer to this transformation simply as DCT. 

Before defining and examining the ACT, some fundamental tools are necessary. In what follows, we introduce the generalized 
Mobius inversion formula tailored for finite series. Additionally, some useful identities for trigonometric sums are also presented. 
We state these preambulary results below. 


Theorem 1 (Generalized Mobius Inversion Formula for Finite Series) Let {/„} be a sequence (e.g., signal samples) such that it 
is nonnull for 1 < n < N and null for n > N. Admit another sequence {g n } defined as 

[ N / n \ 

gn = tik fkn , 

k= 1 

where {a n } is a sequence of scalar coefficients and |_-J is the floor function. Then, 

[ N / n \ 

fn = bmgmn • 

m= 1 

where {b„} is the Dirichlet inverse sequence of{a n }, given that it exists H12V . 
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Proof: It follows by applying the arguments described in the proof of Theorem 3 from 0p. 469] into the Mobius inversion formula 
as shown in |fl3l p. 5561. □ 

If we consider the unitary sequence a n = 1, for n = 1,2,3, ..then the associated Dirichlet inverse sequence is the Mobius sequence 
b„ = n{n), for n = 1,2,3, ... The Mobius function gt(n) fl4i is defined over the positive integers and is given by 


{ 1, if n = 1, 

(—l) 9 , if n can be factorized into q distinct primes, 

0, if n is divisible by a square number. 

In this particular case, arithmetic transform literature simply refers to the above theorem as the Mobius inversion formula for finite 
series Q4). 

The following lemmas are pivotal for subsequent developments. They link usual trigonometric functions to number theoretic 
behavior functions, a connection that is not always trivial mm. 


Lemma 1 (Reed’s Lemma) Let k > 0 and k' >0 be positive integers. Then, 


V („ n I k, ifk\k ', 

2^ cos I 2m—71 I = 


m =0 


0, otherwise, 


and 


*-l / ]f 

^ sin I 2m—K ) = 0. 

m= 0 V k 


We amplified Lemma[I]with the next two results. 

Lemma 2 Let k > 0 and k' > 0 he positive integers. Then, 


V ( k'\ 1 2k, if2k\kf, 

L cos I nm— j = < 

m= o V k J 0, otherwise. 


Proposition 1 Let k > 0 and k' >0 be integers and a be a real quantity. Then, 


k, if k! = 0, 

*-i / jf \ I 

^ cos (m+ a)J = < kcos(27ija), if k\k', k' ± 0, 

0, otherwise. 


Proof: It follows from usual trigonometric manipulations and an application of Lemma\T\ 


□ 


2.2 Arithmetic averages 

Based on existing definitions and concepts inherited from the AFT theory, our initial attempt to define an arithmetic transform 
procedure for the DCT is detailed below. Thus, we consider the following definition for AFT-like averages. 

Definition 1 (Ath AFT-like average) Let the kth average be defined as 

j 2k—i 

*=‘- 2 . 
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Except for the 1 /2-shift, which is naturally imposed by the DCT kernel, the above concept was previously described in OIH UTTI 
as a framework for the AFT. 

Let us investigate the consequences of this definition. Relaxing the integer index constraint of the DCT formulae and substituting 
the inverse DCT relation into the above described A'th average, we must have 


2k—l 


Sk = 


2k 


V„N_ 1 


m =0 

2k— 1 


m =0 


N -1 


= 2k ^ V v £ ak ' Vk ' 


r COS 




A '=0 


2 J 2k—l / N- 1 

E y y o + E COS 1 nm 

m =0 V k f =0 


N 2k 


yv 0 


2 1 vVv ( k ' 

K~ k E E V k' cos ( 

JV Z,C m =0 A '=0 


k = 1,2,... ,1V— 1, 


where 7=1/ \/2 — 1. By interchanging the summations, it holds that 


S k = 


yv 0 


2 1 
N2k 


N -1 2A-1 

L v k' E 

k '=0 m=0 


k' 


cos nm 


Invoking Lemma[2l we have 


k = 1,2,... ,A — 1. 


s ‘ = t/«rVo+,/^I% 


2 1 


JV-l 


7^o- 


A2A: 
2 


^=0 
N-l 

E v/' ■ 

k'=0 


j 2k, if2k\k', 

) 0, otherwise, 

1, if 2£|£', 1 

0, otherwise, J 


A: = 1,2,... ,N— 1. 


Let k 1 = 2sk. Then, 


( 2 ) 


[J fT L^J 

Sk = \ M^ Vo + \h T, y 2 sk, k=l,2,...,N-l. 
v /V v tv s=0 

Without any loss of generality, let us assume that input signals have a null mean value. Therefore, it follows that Vo is zero. This 
assumption does not affect the value of the remaining components of the DCT spectrum. Thus, we find that 


S k 


Va £ V2sk 

v iv J=1 



if v is odd, 
otherwise. 



k=l,2,...,N- 1. 


In order to invert this last expression, we must obtain the Dirichlet inverse of the sequence {0,1,0,1,0,1,...}, as required by 
TheoremQ] However, a sequence can only admit the Dirichlet inverse if and only if its first element is nonnull [ f2j p. 18]. This is 
clearly not the case for the particular sequence under examination. 

We conclude that the derivation of an arithmetic method for the DCT requires more than simply reusing the AFT averages. For 
the DCT computation, specifically tailored averages should be considered. In the next section, such specific derivation is sought. 
However, the mathematical manipulations and arguments shown above prove to be a useful framework. We show next that the 
proposition of DCT specific averages can be greatly benefited from the above exposition. 


2.3 ACT AVERAGES 

In order to derive the ACT, we adjusted the AFT averages as suggested in the following definition. 
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Definition 2 (Ath ACT average) Let the kth ACT average be defined as 



1 t. 1 


where j3 is a fixed real number. 

Considering algebraic manipulations similar to the ones that has led us from Definition [I] to equation (0 and admitting that the 
input signal has null mean, mutatis mutandis, we obtain 



S k 


Thus, invoking Proposition!]] we have 



Performing the substitution k! = ks, it follows that 



Yj cos(2jtsP)V sk , k= 1,2,...,A- 1. 


S k 


Observe that the above expression is suitable for the application of the generalized Mobius inversion formula for finite series. Under 
the notation of TheoremQ] we recognize a n = cos(27Tnj3), n = 1,2,3,... 

Incidentally, not always the Dirichlet inverse of { a n } is well-defined. Only when a\ / 0, the existence of the Dirichlet inverse 
can be considered El- Thus, we must impose cos(27 z[5) / 0 as a necessary condition for the derivation of the ACT procedure. 
This issue is precisely the point that Definition[I]misses. 

However, finding the Dirichlet inverse of {«„}, say {b„ } for arbitrary values of [5 can be a cumbersome maneuver. Therefore, 
we separated two particular useful cases: (i) /3 = 0 and (ii) /3 = 1 /2. Notice that /3 = 1/4 leads to a non-invertible sequence, since 
this makes a i = 0. 

For [5 = 0, we have the unitary sequence a n = 1 and b„ = jJ,(n), for n = 1,2,3,... This is usually the situation addressed in 
standard AFT analysis. On the other hand, setting /3 = 1/2 yields a n = (—1)", n = 1,2,3,... In this case, the Dirichlet inverse is 
not immediately recognized, but is can be obtained analytically. In the Appendix, we derive the sought Dirichlet inverse, which is 
given by 



— fl(n ), if n is odd, 

—2 m-1 jti(2 ~ m n), if n = 2 m s, where s is odd. 


(3) 


The first 32 terms of {b,,} are listed below: 


-1, -1, 1, -2, 1, 1, 1, -4, 

0 , 1 , 1 , 2 , 1 , 1 , - 1 , - 8 , 

1, 0, 1, 2, -1, 1, 1, 4, 

0, 1, 0, 2, 1, -1, 1, —16. 


In the context of digital signal processing this is a potentially useful sequence, since multiplying a given number by a power-of-two 
can be implemented by simple bit shift operations, which possess a low computational complexity. 
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Regardless the choice of /l, to invert the following expression given by 


S k = \l^ £ a s V sk , k= 1,2,...,1V— 1, 
v 7V j=i 

we may consider an auxiliary sequence given by G k = \J2/NV k for all k. Thus, a direct application of the generalized Mobius 
inversion formula for finite series just tells us that 


L^t 1 ] 

Gk= £ b ‘ S kh k= 1,2,... ,N — l, 

i=i 

where {£>„} is the Dirichlet inverse of {a,,}. Finally, undoing the auxiliary substitution we can write 

Vk = \ ^r t b ‘ S kh k = 1,2,1. (4) 

v 1 1=1 

Depending on whether {a „} is the unitary sequence or the alternate sequence —1,1, —1,1,..., the inverse sequence {/?„} is the 
Mobius sequence or the sequence described in (0), respectively. 

Recognizing that unitary and Mobius sequences constitute a simpler pair of Dirichlet inverse sequences, in the rest of this paper 
we focus our attention to the case P = 0. However, all ensuing developments encompass the proposed alternative formulation 
(j3 = 1 /2) without significant modifications. Nevertheless, the case p = 1 /2 can furnish a link between Definitions Q] and 0 Indeed, 
for P = 1/2, we have 


1 k ~ l 
S k = ~ k L v 

K m =0 


2(m+'' w 1 




1 k ~ X 

k ^ ' ; (2m+l) 
K m=0 


N 1 


1 2k-1 


m= 1 
m odd 


This last summation corresponds to the average shown in Definition0 when only odd values of the dummy index are considered. 
In other words, if the samples required by Definition0are submitted to a downsampling by a factor of 2, then Definition0collapses 
to Definition0for p = 1/2. In a sense, this observation establishes a connection between both formalisms. 

As far as the computational complexity of the Mobius inversion formulae are concerned, we can provide the following probabilis¬ 
tic reasoning. The probability that a randomly chosen integer is not divisible by a perfect square is 6/k 2 « 0.61 lfl2l p. 4]. Therefore, 
61% of the values of the Mobius function are zeros; meaning that the computation of G k (or V k ) requires (1 ■ 6/n 2 ) (A 1) /k\ 
additions/subtractions on average. 

Up to this point, we assumed that input signals have null mean. Now, let us remove this restriction. Let v = ^ L^=o v « the 
mean value of the input signal. Thus, an arbitrary signal v' can be converted into a null mean signal by simply subtracting the 
quantity v 7 . Therefore, the kth ACT averages S' k associated to V — v can be manipulated as follows: 



k—l 


E(v 

ra=0 


2mf 


1 

1 


V) 




= S k -v, k = 1,2,... ,N — l. 
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Consequently, equation ([4} can be rearranged, yielding 


V k 





Considering the Mertens function M(n), which is defined as M(ri) = YJln=\ jJ.(m) lfj~8l p. 272], we obtain a more compact 
expression as 


WJ 


Vk=\l 2 ( I 


N- 1 


, k = 1,2,... ,N — 


(5) 


The zeroth component Vq is expressed separately by Vq = y/Nv. 


3 ACT INTERPOLATION 
3.1 Fractional indices manipulation 

Now we must address the problem of handling fractional index samples. Usually, arithmetic transform literature addresses this issue 
by prescribing the utilization of zero- or first-order interpolation methods. We argue that such a naive approach is not the proper 
way of interpolating. 

Again, let us relax the integer index constraint of the DCT formulae. Considering a possibly noninteger quantity r , the following 
construction is derived given by 

T V (nk{r + 1/2) 

Vr = \/ 77 L «^ cos ' 


k =0 


N 


Taking into account the direct transformation formula, we obtain 


W /TV (nk(n + ll 2)\\ (nk(r+\/2) 

Vr = \ — a k\ \ T-jUk Vn COS | -—- 1 I COS ' 

\ V A n=0 


N 


N 


Inverting the summation order yields 


N n =0 k =0 

Then, let us define the ACT weighting function as 


2 V' V 2 fnk{n+ 1/2) \ (nk{r+ 1/2) 

Vr = — L Vn L a k COS I -- I cos 


N 


N 


Mr) 4 1 £ eel™ (aOa+im - <■ ” l(r+1/2) 


k=0 

1 2 


N 


N 


N—l 


« ^ 


COS 


nk(n + 1/2) \ /^(r+1/2) 


N N£ 0 \ N J V A 


cos 


n = 0,1,... ,A— 1. 


Thus, the samples associated to fractional indices can be obtained after a linear combination of the available uniformly obtained 
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samples. Hence, 


N-l 

V r = E w n(r)v n . 

n =0 


Notice also that if r is integer, we have that 


w n (r) 


1, if n = r, 

0, otherwise. 


( 6 ) 


This fact stems from the orthogonality properties of the transformation kernel. 

We further investigate the suggested weighting function. The next proposition states that the discussed weighting functions are 
indeed inherently normalized. 


Proposition 2 The ACT weighting function satisfies the following summation formula 

N-l 

E w «( r ) = L 

72=0 


Proof: Usual trigonometric manipulations furnish the derivations below: 

'Jtk(n+ 1/2) 


N-l N-l / i 9 N-l 

E w «W = E (-^+^ E cos 


n =0 


77=0 

= -i + 


N 


k =0 

2 N-l / I 1 /T\ \ N-l 


COS 


nk(r + 1/2) 


N 


cos 


k=0 


nk(r+ 1/2) 
N 


E cos 

72=0 


N 

nk{n+ 1/2) 
N 


However, for each k, the inner summation can be expanded as follows: 


N-l 

E cos 


77=0 


7 ik{n + 1 /2) 
N 


nk 
2 N 


N-l 


E cos ( 

n =0 V 


/ 2 n(k/ 2 )n 


N 


nk 

2N 


N-l 

E sin 

72=0 


. / 2 n(k/ 2 )n 


V AT 


Since dummy index k runs from 0 to N — 1, we have that N never divides k/2. Therefore, applying Lemma\T\ we obtain that 

nk(n + 1/2)' 


N-l 

E cos 

72=0 


N 


( nk\ N, ifk = 0, 

= cos - < 

\2NJ 0, otherwise, 

_ [ N, ifk = 0, 

0, otherwise. 


Finally, returning to the previous double summation, we establish that 


V ,2V (nk{r+\/2)\) N, ifk = 0, 

E Mr) = -l + jj E cos ' 


72=0 


k =0 


N 


0, otherwise. 


= 1. 


□ 

Regardless of the considered block-length, the use of the ACT interpolation results in an exact calculation of the DCT spectrum. 
Indeed, no approximation was considered in any of our arguments. 
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Figure 1: Block diagram of the 8-point ACT. Dashed lines indicate multiplication by —1. 

3.2 An example: the 8-point ACT 

To illustrate the ACT structure and its interpolation scheme, we devised an example. Intentionally, we selected the short block- 
length transformation furnished by the 8-point DCT. This particular block-length is widely adopted in several image and video 
coding standards, such as JPEG, MPEG-1, MPEG-2, H.261, and H.263 |fl9l . The 8-point DCT is also subject to an extensive 
analysis in HD. In the following, we set j3 = 0. 

The first step of the ACT procedure consists of the identification of the necessary interpolation points. According to Definition^ 
these points are given by 2m | — ^ for k = 1,2,... ,7 and m = 0,1,2,... k— 1. Therefore, we find the following fractional indices: r £ 
{ — 2,l!>T>Tn,2,IJ>T’fo’fl>T}- P'S- Qjdepicts a block diagram of the full algorithm. The ACT interpolation block calculates 
the required samples v r according to the discussed interpolation procedure. 

The fractional index samples are then employed to obtain the ACT averages. The block of ACT averages simply implements 
the following set of equations: 


*Si = v i 
2 

25*2 = Si + Vl5 

3 S 3 = S\ 2v 29 
TO 

4 S 4 = S 2 2v 7 

5 S 5 = Si + 2v 27 + 2v 59 
TO TO 

6 S 6 = S 2 H“ S 3 + 2v 13 
TO 

1S~] = +2V25 + 2v57 + 2v89 . 

U I? U 

Finally, the ACT averages are combined with respect to the Mobius function (cf. Qi). The resulting calculation involves no approx¬ 
imations and furnishes the exact DCT spectrum. 


3.3 Asymptotic analysis 

In this section, we closely examine the weighting function required by the interpolation procedure. We aim to propose simpler 
interpolation expressions allowing an efficient computation of the ACT. 

First, we note that the ACT weighting function can be formulated in closed form as detailed below. In fact, invoking elementary 
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trigonometric identities, we can establish the following relations: 


N -1 


v„(r) = —i— y 


A? A 

1 1 
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k=0 
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/t=0 


COS 


1 l^ 1 

-1-> cos 

N A 
Ly ly k =0 


nk(n +1 / 2 ) 

TV 

7tA:(« + r+ 1) 
TV 

n k(n + r+ l] 

TV 


cos 


nk(r+1 /2) 

TV 

7tA:(h — r) 


-cos 


i JV-l 

+ — y cos 
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k=0 


A 

nk{n — r) 
TV 


, « = 0,1,... ,TV — 1. 


The above trigonometric summations can be given in terms of the Dirichlet kernel [20] p. 312], Therefore, it holds that 


w„(r) = - — + 


1 1/11 


TV TV \2 2 

2TV 


( n < 

— In 

\N y 


+ l ( — (w + r + 1 


+ 


1 /I 1 


TV \ 2 


■ + 7%- 1 


(s ( "- r) ) 


f (oA.- | (i ( “ + r+l ) )+D„_,(^(“- r ) )) 


where D N (x) = Sln ^^^' x) denotes the Dirichlet kernel. 

The similarity between the Dirichlet kernel and the sine function is apparent. Indeed, as TV —> °°, these functions are interchange¬ 
able in several situations GQ- However, even for small values of TV such an approximation is good GQ p. 180]. For instance, 
taking A = 8 and centered functions at the origin, it follows that the mean square error (MSE) of implied approximation is less than 
2 x 1 0 3 over the interval [— TV/2, TV/2]. Thus, the discussed weighting function assumes a limiting form given by the sum of two 
translated sampling functions as 


lim w n (r) = sinc(n + r + 1) + sinc(n — r), 

N-> °° 


. / \ A smf.T.v) 

where sinc(x) = —^- L . 

Notice that this asymptotic expression connects the ACT interpolation to the sine function interpolation. Signal interpolation 
according to the sine function can be efficiently implemented in the time domain 112211231 . Additionally, fractional delay FIR filtering 
methods offer another computational approach 112411251 . 

The limiting form of w n {r) leads us to draw some additional conclusions on its asymptotic behavior. Let [•] denote the nearest 
integer function as implemented in C or Matlab programming languages and r be a fractional interpolation point. We observed that, 
when 0 < [r] < TV — 1, the asymptotic weighting function is essentially governed by a single sine function. The role of sine (n + r + 1) 
could be neglect, since its argument would be large enough. Indeed, this function approaches to zero according to (9(1 /n). Thus, in 
this case we say that 


lim w„(r) « sinc(« — r ). 

iV->°o 

Fig.[2]shows some plots of sine (/i — r) for several values of TV. These plots intuitively indicate the use of a linear approximation 
for the main lobe of the sine function. Additionally, we assumed that the effect of the remaining lobes is negligible. We also observed 
that, for [r] £ {— 1, TV}, both sine functions are relevant since they overlap each other. This last situation was treated separately. 

In terms of the above discussion, we derived an empirical approximation procedure that considers at most two uniform samples 
to render an interpolated sample. In Fig. [3] algorithmic details of the suggested approximate interpolation are given. In the proposed 
heuristic algorithm, we admit an auxiliary quantity A = r— \r], an error tolerance e, and a scaling factor a ft: 1.2, for A = 8. This 
last quantity was found according to a standard linear fitting procedure. 

Considering a conservative choice of e = 0.1, we employed the proposed heuristics to calculate the approximate DCT spectra of 
256 randomly generated signal vectors. A short block-length A = 8 was deliberately selected. The elements of the input signal were 
chosen to be distributed according to a standard uniform distribution. The resulting average MSE due to the discussed approximation 
was as low as 4.7 x 10 / This figure is comparable to the MSE associated to some integer approximation algorithms for the 8-point 
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Input: Fractional number r, tolerance £, block-length size N. 
Output: Approximate values of the ACT weighting function. 
Method: Heuristic approximation. 


a <— 1.2, A <— r— [r], w«— 0 

if |A| < £ then 

W [r] = 1 

return 
end if 

if [r] € [1, JV — 2] then 

W[r]-l <- (l A l — A )/ 2 

W[ r ] <- 1 - ! A I 

w [r]+l (l A l + A )/ 2 

else if [r] = 0 then 

wu «— 1 — |A| 

M '[r] + 1 (l A l + A )/2 

else if [r] = N — 1 then 

M'H-1 <- (l A l — A )/2 

W{r] 1 - l A l 

else if [r] = — 1 then 

wo 1 

wi <-0.35 

else if [r] = N then 

wjv-2 ^-0.35 

WjV —i ^— 1 

end if 

return a ■ w 


Figure 3: Algorithm for computing approximate values of the weighting function w, (r), i = 0,1,... ,N — 1. 

DCT (e.g., integer cosine transform, C-matrix transform) as detailed in HU- 

4 Matrix-vector representation 

Further insight on the nature of the ACT can be obtained when the previous constructions are represented in matrix-vector 

lT 


form. Let the input signal and its associated DCT spectrum be denoted by column vectors v = 

1 T 


I'O, Vl, 


Lv-1 


and 


V = 


Vo, Vi, 


Vv- 


, respectively. Additionally, consider the DCT matrix C, whose elements are defined accord¬ 


ing to <□}: = -\/2/N(XkCos(nk{n + 1/2 )/N), for k,n = 0,1 1. The above quantities are related by V = C • v and 

C^ 1 =C T QUp.41], 

In order to render the ACT structure in matrix-vector form, we need to introduce some special matrices. 

Definition 3 (Mobius matrix) The ( i,j) element of the N-order Mobius matrix M,.y is given by Li ( j/i), whenever i divides j, for 
i,j= 1,2,..., N. Otherwise, it is zero. 

By construction, the Mobius matrix is upper triangular with unity diagonal elements. Thus, it is always non-singular and its 
determinant is unity for any dimension. The inverse of the Mobius matrix can be directly obtained without calling an inversion 
procedure. The (i.j) element of M 1 is 1, if i divides /; otherwise, it is zero. This fact stems from the Mobius inversion formula 
relations. 

Additionally, we consider the extended Mobius matrix defined by 


M' = diag(l,M^_ 1 ) = 


1 0 
0 


0 


M, 


JV—t 


where the operator diag(-) returns a diagonal matrix. 
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We may also admit a vector of averages expressed by S = 


So, Si, 


Sn- 


where the zeroth average is defined 


separately as So — v i’ anc * t ^ le rema i n ing values are the Ath ACT averages for k = 1,2..... /V 1. Notice that DC component 

of the spectrum is related to the zeroth average Vo = y 7 N/2Sq. 

The ACT framework can provide a new insight into the DCT spectrum. In fact, equation (0 indicates that the DCT spectrum 
can be separated into two parts: (i) one due to the Mobius combination of the ACT averages and (ii) another due to the Mertens 
function. Let these two parts be termed the Mobius and the Mertens parts of the DCT spectrum denoted by Vi and V 2 , respectively. 
Thus, the DCT spectrum can be decomposed as V = Vi + V 2 . 

Joining the above structures it follows that the vector S is related to the Mobius part of the DCT spectrum via the Mobius 
extended matrix according to 


Vi = 


•M'-S. 


(7) 


In the light of the proposed ACT interpolation, we can recast the ACT averages in terms of the weighting function. Thus, 
invoking 0, we have 


1 k ~ 1 

1 k ~ ]1 ( n - 1 ( N 1 

= ^L L v (1 vv„ 2m--- 
K m =0 \n=0 V K Z 

Inverting the order of the summations, we may also obtain 




N-l 


1 


k -1 


Sk = L 7 L w "( 2m T - 


n =0 


ra=0 


N 1 




For a fixed n, the inner expression in parenthesis is an average of particular weighting values at different interpolating points. This 
term depends only on n and k being independent of the input vector v. Then, we can separate this expression for better understanding. 
Let us define the Ath weighting average as 


A it, 1 / A 1 

W kn = 7 I I 2m— - - 


m =0 


k=l,...,N-l, n = 0,...,N-l. 


These averages give rise to the construction of the following matrix: 


W =\W k , n ], k= — 1,« = 0,...,A — 1. 

Since the rows of W are convex combinations of the weighting function, Proposition[2]indicates that the sum of elements across 
rows of W is equal to 1, i.e., Y,n=o W k ,n = 1, for A = 1..... /V — 1. By augmenting the matrix W with the inclusion of an extra row 
and padding, as described below. 


we may write the following relation: 


1 /N 1 /N ■■■ l/N 
0 

w 


0 


S = a • W' • v, 
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where a = diag(ao,ai,...i). 

An application of this last expression into 0 furnishes the final matrix-vector form for the Mobius part of the considered 
procedure. Hence, 


Vi 



•M'-a-W'-v. 


Effectively, the transformation matrix that relates Vi and v is given by 


Ci = \/ — - M' • a-W'. 


Now considering the spectral part V 2 due to the Mertens function in 0, we advance the following additional matrix: 


C 2 




N— 1 
1 


M 


N- 1 
2 


.,M 


N- 1 
N- 1 


• Ijv 1 


T 

N’ 


where 1# is a column vector consisting of unit elements. Of course, if the input signal has null mean, then the Mertens part of the 
DCT spectrum is always zero: C 2 • v = 0. 

In view of the above, we must have 


V = Vi + v 2 
= Cl -v + C 2 v 
= (Ci + C 2 )-v. 

This manipulation suggests that the DCT matrix C can be decomposed as C = Ci + C 2 . Notice that when the input signal has null 
mean, a new formulation for the DCT matrix arises. In this case, we can establish that 

C v = Ci • v. 

This alternative transformation matrix for the DCT can be considered as a starting point to derive new algorithms under the constraint 
of null mean input signals. 


5 Conclusion and Remarks 

The main contributions of this paper are (i) a new arithmetic transform and (ii) the elucidation of the arithmetic transform interpo¬ 
lation issues. The introduced arithmetic cosine transform is a number-theoretic based algorithm devoted to the DCT computation. 
Therefore, DSP applications that require a DCT evaluation are potential candidates for the use of the ACT method. 

Differently from the standard AFT analysis, we generalized existing inversion formula, allowing new frameworks for the arith¬ 
metic transform theory. Besides the Mobius sequence, we identified alternative adequate sequences for the suggested method. 

Moreover, the introduced interpolation method allows the exact computation of the DCT spectrum, even for small block-lengths. 
This is particularly distinct from the existing AFT algorithms, which (i) inevitably introduce approximation errors and (ii) tend 
to excel only when large block-lengths are considered. Using arithmetic methods (e.g., AFT) for small block-length transform 
evaluation was a challenging issue, for which area literature could not furnish adequate solutions until now. 

The complexity of the ACT procedure is mainly due to the interpolation step. If the sampling process could be adjusted in 
order to collect samples natively in an appropriate nonuniform fashion, then the ACT computation would not need any sort of 
interpolation. Indeed, the DCT would be exactly computed after some few additions/subtractions associated to the Mobius function. 
Other than that, the arithmetic transform philosophy can concentrate the transform computational complexity into the interpolation 
block. This is a different perspective for the design of transform procedures. For instance, fast algorithms for fractional delay 
filtering Ii24j now receive an additional motivation under the arithmetic cosine transform paradigm. 
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Finally, the existence of efficient frameworks for bidimensional AFTs suggests a venue for extending the proposed ACT into 
the bidimensional case |10||26fj28l . Without much effort, the proposed method could be used as the fundamental building block for 
a bidimensional arithmetic cosine transform. On the other hand, a dedicate analysis for the bidimensional case could prove to be 
more appropriate. In particular, an application Feig-Winograd direct approach could avoid row-column methods ||29l . This could be 
a more attractive method for the bidimensional ACT. In any case, the characterization of the exact interpolation for bidimensional 
signals is an open topic. 
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A The Dirichlet inverse of {(-1)"} 

In order to derive the Dirichlet inverse of the sequence a„ = (— 1)", for n = 1.2.3,..., let us examine its associated Dirichlet series. 
A result from the theory of functions [ (301 p. 337] states that 

“ f—11" 

L L -r-=-(l-2 1 " , )CW, 9tM>0, 

n= 1 n 

where £( 5 ) = £“ =1 1 /n s is the Riemann zeta function. Therefore, we directly have that the closed form of Dirichlet series of {«„}, 
A{s), is 


A ( S) = _( 1 _ 2 1 -)C(,). 


The Dirichlet inverse of {a n } is a sequence {b n } such that its Dirichlet series, B{s) = LJT=i b n /n s , is equal to 1 /A(s). Thus, we can 
maintain that 


B{s) = 


1 1 
l-2'- s CW 


1 y> /j( ? 0 

1 — 2 l ~ s n s 

n = 1 


Before finding {b n }, we must identify the Dirichlet series of 1/(1 — 2 1 '). This is necessary in order to put B(s) as a product of two 
Dirichlet series and apply the convolution theorem for Dirichlet series. Accordingly, we have that 


1 

1 - 2 1 -* 


E( 2 l ~ s ) k = 


k =0 


OO 


E 

A -=0 


2 k 

(2 k y ' 
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This final expression is already in the Dirichlet series format. Therefore, the sequence associated to 1/(1 — 2 1 s ) is simply 


Cn — 


n, if n is a power of two, 
0 , otherwise. 


Returning to the expression for B(s), we can write as 


bw=- £ £ " w 


I 

n= 1 


(c®n)(n) 


where © denotes the Dirichlet convolution. By the equivalence property of Dirichlet series, we conclude that b n = —( c® ). 

Now let us evaluate ( c®ji)(n). Observe that if n is odd, then the only divisor of n which is a power of two is the unit. Therefore, 
we obtain 


(c® n){n) = Y, c dH(n/d) = cin(n) = n(n). 

d\n 


On the other hand, admit that n is even in the form n = 2"',v, where s is an odd integer and m is a positive integer. Considering that 
(i) the Dirichlet convolution of two multiplicative functions results in a multiplicative function dp. 35] and (ii) the sequence {c„} 
is multiplicative with respect to n (a fairly direct result), we maintain that 

(< c®n)(n ) = (c®n)(2 m ) ■ ( c®n)(s ) = (c©ju)(2 m ) -ju(s). 

According to the definition of the Dirichlet convolution, the expansion of (c® jl)(2 m ) yields 

(c® /J,)(2 m ) = cifi(2 m ) T C 2 B(2 m *) 4-f c 2 m-i/t(2) +C 2 m /i(l). 


Due to the Mobius function, only the last two terms are possibly nonnull. Thus, we have 

(c ©;U)(2 m ) = c 2m _, At (2) + c 2 mH (1) = 2 m ~ 1 (-1) + 2 m ( 1) = 2 m ~ 1 . 

Joining the above manipulations, we have that 


(■ c®n)(n ) 


M(«), 

2 m ~ 1 n(2~ m n), 


if n is odd, 

if n = 2 m s, where s is odd. 
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